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Locating the global minimum of a complex potential en- 
ergy surface is facilitated by considering a homotopy, namely 
a family of surfaces that interpolate continuously from an ar- 
bitrary initial potential to the system under consideration. 
Different strategies can be used to follow the evolving min- 
ima. It is possible to enhance the probability of locating the 
global minimum through a heuristic choice of interpolation 
schemes and parameters, and the continuously evolving po- 
tential landscape reduces the probability of trapping in lo- 
cal minima. In application to a model problem, finding the 
ground-state configuration and energy of rare-gas (Lennard- 
Jones) atomic clusters, we demonstrate the utility and efficacy 
of this method. 



Introduction: Global optimization problems can of- 
ten be formulated in terms of finding the minimum (or 
maximum) of a multidimensional potential energy sur- 
face (PES) . Such problems, which occur in a variety of ar- 
eas, are of considerable practical and theoretical interest 
1^ . The "energy landscape" |^ paradigm is particularly 
useful when the potential energy function is continuously 
varying with the physical configurations relevant to the 
problem . An example of such a situation is the protein- 
folding problem namely determining the native con- 
figuration of complex molecules given their atomic com- 
position. A simpler variant is the determination of the 
ground state configuration of atomic or molecular clus- 
ters |l. 

In this Letter we propose a new homotopy method to 
study such problems by a controlled deformation of the 
potential energy surface. If Vf is the potential energy hy- 
persurface under consideration, we study the landscapes 



V(a) = {l~a)Vi+aVf, 



(1) 



with a a parameter. Given a choice of initial potential, 
Vi, this is a 1-parameter family of potential energy sur- 
faces which smoothly evolves from Vi into Vf as a varies 
from — > 1. 

The minima of the landscapes continuously change 
with a, and in order to track them, one of two strate- 
gies are possible. Varying the interpolation parameter a 
in a finite number of steps, a standard technique such 
as conjugate gradient (CG) minimization Q can be em- 
ployed at each a. On the other hand, one can consider a 
as a time-dependent function such that the PES evolves 
according to 



V{t) = {l-h{t))V, + h{t)Vf, 



(2) 



where h{t) is suitably chosen with h{Q) ~ 0, and 
\\mt->Th{t) 1. Over a timescale T, therefore, the po- 
tential deforms from the initial to the desired potential 
energy surface, and the evolving minima can be tracked, 
for example, by following the damped dynamics in this 
potential via molecular dynamics (MD) simulation. 

In the present work we follow both these strategies, 
and show how homotopic deformation facilitates loca- 
tion of the global minimum in a model problem. Sim- 
ilar (so-called "continuation") homotopic methods have 
frequently been employed in related situations, as for ex- 
ample in finding roots of polynomial equations in several 
variables |Q or in the mean-field dynamics in attractor 
neural-networks 

Different global optimization methods frequently find 
optimal solutions by elimination, by seeking lower and 
lower minima. Trapping in local minima — and escape 
from these minima — is a major practical issue. A num- 
ber of different strategies have been suggested in order to 
engineer escape from local minima. These include both 
techniques to allow for large excursions in the phase space 
by the use of temperature or similar auxiliary parameters 
(such as simulated annealing [|l0| and its variants ) 
as well as methods that deform the potential energy sur- 
face. The diffusion equation method ||l^ and the distance 
scaling method |Q fall in this latter class. Other meth- 
ods utilize both strategies, as for example the stochastic 
tunneling method where simulated annealing is per- 
formed on a surface where the barriers are exponentially 
reduced so as to facilitate escape from local minima, the 
landscape paying technique p6| , or the basin hopping 
technique which replaces the potential surface by a 
set of piecewise fiat regions. 

The present technique is in the class of optimization 
methods that exploit potential surface deformation to 
avoid trapping in local minima. The interpolation pa- 
rameter a, or the switching functions h{t) smoothly con- 
vert one PES into another. The intermediate potentials 
are qualitatively not very different from the asymptotic 
potential in terms of the number of minima and maxima, 
although the relative depths and curvatures are quite dif- 
ferent. As we discuss below, this feature contributes to 
efficiency of the present technique in locating minima. 
The lowest energy achieved when an ensemble of suit- 
ably compact initial configurations is evolved is taken as 
the ground state prediction of this method. 

Application: The problem of minimum energy con- 
figuration determination for N particle atomic clusters 
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is computationally hard, and the validity of a global so- 
lution cannot, typically, be verified. Existing data for 
global minima are usually the "lowest minima as yet 
located" in all but the simplest cases. A variety of global 
optimization techniques have been applied to this prob- 
lem [ p7|Jl8| with differing degrees of success. 

For the most extensively studied such systems, namely 
model rare-gas clusters, the potential energy surface 
(PES) is an additive pairwise Lennard-Jones interaction, 

where Vij is the distance between particles i and j, and 
e, a are the standard Lennard-Jones parameters. The po- 
tential energy landscape varies greatly with cluster size. 
Notable difficult optimization problems in this regard 
are, for example, 38, 75, or 98 atom clusters, where the 
potential energy surface has the so-called multiple funnel 
structure (20[|2l|,||. 

In the implementation of the MD approach we proceed 
as follows. Vi is taken to be a pairwise sum of harmonic 
terms V{rij) — {vij — 2^/^cr)^/2. We perform molecular 
dynamics simulations [^3| of the N particle system, with 
an additional damping term for each particle, 

m?^+"^^^ + ^^ ^ 0,1^1,..., N (4) 

where m is the mass of the particle and 7 is the damping 
coefficient. The internal timescale of interparticle vibra- 
tions depends on the parameters m, a and e. For a given 
switching function h{t) (we have explored a variety of 
such functions listed in Table I) the adiabatic time scale 
is set by the parameter C; the entire system dynamics 
thus has two external time scales C,~^ and m7~^. In the 
limit 7 — > 00, our procedure reduces to a steepest descent 
minimization on the evolving potential. The dynamics of 
the system is followed until a stationary configuration is 
reached. 

In order to quantitatively assess the efficiency of this 
procedure, we define the measure 

^ Number of ground state configurations 
^ Total number of condensates 

a condensate being a configuration such that all atoms 
are within a single cluster. This is clearly a function of 7 
and C,. For the ground state energy, comparison is made 
to the existing benchmark calculations already available 
for Lennard-Jones clusters [ p^ . 

In the CG approach, Vi is taken to be /3X)jLi(^j ~^)^' 
being the (random) initial position for the jth atom. 
This choice of Vi ensures that the initial configuration is 
the exact global minimum for the potential energy sur- 
face, Eq. (|^) with a — Q; P \s a. constant that tunes 
the curvature of the PES. The parameter a is then var- 
ied from to 1 in Ns discrete steps; the result of the 
CG minimization (we follow the standard method [|| ) at 



each step is taken to be the starting configuration for the 
CG minimization at the next value of a. In this latter 
approach, therefore, the attempt is to allow the global 
minimum itself to evolve homotopically. 

Results: The present application is intended to be illus- 
trative rather than exhaustive. We have systematically 
studied different cluster sizes up to TV = 40 and in all 
cases the calculated ground-state energy and configura- 
tion matches existing results exactly. This includes the 
difficult case of the 38-atom cluster which is an interest- 
ing and important test for any optimization method [ pO| . 
The number of minima increases exponentially with clus- 
ter size p^ ; for LJ7 there are 4 minima, while for LJ55 
the number exceeds 10^° [^. Detailed results, which clar- 
ify some aspects of the present technique are presented 
for the cases of iV=19,22 (MD) and A^=38 (CG). 

In the MD version of the present technique, in the ab- 
sence of switching, namely in the sudden limit V{t) — Vf, 
the system quickly settles into the nearest available mini- 
mum based on the level of damping introduced. By start- 
ing from an ensemble of initial conditions, a variety of dif- 
ferent minima are reached but the probability of finding 
the true ground state is essentially zero for large clusters. 
With an adiabatic switch , the results are dramatically 
different. The continuous evolution of the potential en- 
ergy landscape is a key factor in permitting escape from 
local minima. Only asymptotically does the system come 
to rest, but until then, there is always residual kinetic 
energy due to which the system avoids being trapped by 
small barriers. Shown in Fig. 1 is the typical variation of 
potential energy (in units of e), which is nonmonotonic 
once the adiabatic switching is incorporated. Regardless 
of the actual form of the switching, more than 85% of all 
initially random configurations condense, except in the 
case where the switching is applied to the repulsive term 
of the potential. Representative data is given in Table I. 

As emphasized, the adiabatic optimization proposed 
here is heuristic. The optimal choice for the parameters 
7, C for a given cluster size depend on a number of fea- 
tures such as the interaction potential parameters and the 
inherent time-scales. By scanning over reasonable values 
of the parameters, it is possible to determine regions in 
parameter space with a higher than average probability 
of reaching the ground state. It also appears that adi- 
abaticity is crucial since the probability of reaching the 
ground state increases substantially with decreasing C: 
Vg is shown versus C for the 19-atom case in Fig. 2. 

In the CG method of following minima during homo- 
topy, the probability of reaching the global minimum is 
enhanced through the modification of the PES curvature. 
Since Vi adds a uniform positive curvature at the interme- 
diate stages it effectively suppresses or eliminates many 
barrier. To perform some benchmarking of the advan- 
tage this gives, we present, in Table II, data pertaining 
to finding the global minimum for LJ38 comparing the 
present method and the basin-hopping technique. The 
three lowest minima are at energies -173.928, -173.252 
and -173.134 respectively. In either method, all particles 
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are initially placed randomly inside a sphere of radius 
i^y^^- III the parameter (3= 100. In our implemen- 
tation of the basin-hopping algorithm, coordinate dis- 
placements are random in the interval [-0.3,0.3] and the 
temperature is taken to be 2. An overall confining poten- 
tial of the form Vc = exp(20(ri - a)), a = 1 + (f 
was added to prevent dissociation, and a standard Polak- 
Ribiere algorithm was used for the conjugate-gradient 
minimization with tolerance set between 10~^ and 
10~^. The average computational effort required is a 
product of the number of trials needed in order to get to 
the ground state on average and the number of function 
and derivative calls per trial. In our implementation of 
the algorithms, we find that the reduction in computa- 
tional effort in locating the global minimum through the 
homotopy method is about 40%. The relative efficiencies 
can, however, vary depending on the actual choice of the 
various adjustable parameters in the two techniques. In 
either the MD or the CG version, configurations that do 
not reach the global minimum still typically tend to find 
the lowest energy states, so that a by-product of this 
methodology is a considerably detailed map of the low 
excitation regime of the cluster. This feature, however, 
is not unique to the present method. 

Summary: We have presented here a method for global 
optimization which relies on the guided evolution of the 
underlying landscape. The methodology for finding min- 
ima on this surface can vary, and in the examples pre- 
sented here, we have used both the conjugate gradient 
technique as well as damped molecular dynamics. (Dy- 
namics in the landscape has been incorporated in other 
techniques, for example in genetic algorithms ||2^.) As in 
other methods, apart from the global minimum, we also 
obtain a detailed picture of the excitation spectrum. 

Within the context of cluster geometry determination 
itself, several issues need to be addressed. The adia- 
batic method can be shown to locate ground states even 
when there are bifurcations along the deformation path- 
way ||2^ . Is it possible to design more efficient homotopic 
deformations? What is the role of Vi in controlling the 
efficiency? 

The application here, though in some ways a model 
problem, has all the complications that arise in more gen- 
eral optimization problems. The success of this simple 
technique is therefore encouraging. 
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Table I: Representative results using the MD version of 
the homotopy method for LJig and LJ22 as examples of 
magic and nonmagic clusters. The parameters are C = 7 
= 0.5, for 1000 trials starting from random initial condi- 
tions. 



3 





P 


9 


h{t) 


LJig 


LJ22 


1 (No Switching) 


0.0 


0.0 


1 - exp(-Ct) 


0.0081 


0.0083 


sm 


0.0065 


0.0181 


T 


0.0044 


0.0138 


(tanh(C< - 10) + l)/2 


0.0124 


0.0176 


1 -exp(-C^) cos^(3Ct) 


0.0111 


0.0041 



Table II: Comparative analysis of the homotopy 
method and basin-hopping. For each method, initial 
configurations are evolved to find the global minimum in 
100 instances for the LJ38 cluster. Nj,j = 0, 1, 2 are the 
number of times the lowest three minima are found in 
the two methods; the number of function and derivative 
calls needed (per initial condition) are also indicated to 
give an estimate of the computational effort involved. 



Optimization 
Method 


Nr 
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Ni 


N2 


Function 
Calls 


Derivative 
Calls 


Basin 
Hopping 


937674 


100 


239 


941 


3495 


154 


Homotopy 
Method 


195690 


100 


4 


173 


9260 


475 



0.016 



0.012 
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FIG. 2. Probability of reaching the ground state, Vg, as a 
function of C for h{t) = 1 — exp —(^t, for the cluster LJ19. 
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FIG. 1. Typical variation of potential energy (in units of 
e) with time for the condensation of L J22 , for the case of no 
switching, h{t) = 1 (dashed line), and with switching (solid 
line) using h{t) = 1 — exp(— ("f). 
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